2.3. Resonances of a 2d dumbbell

by M. Wess, 2024

This Notebook is part of the td_evp documentation on the implementation of time-domain methods for resonance problems in NGSolve.

Introduction

This notebook is designed to present an implementation of the algorithm for solving resonance problems presented in

• [NW24a] Lothar Nannen and Markus Wess. A krylov eigenvalue solver based on filtered time domain solutions. 2024. arXiv:2402.08515.

In particular, we solve the two dimensional problems from Section 3.1.

Implementation

We start bei doing the necessary imports

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *

import numpy as np
import scipy.linalg as spl
import matplotlib.pyplot as pl

The filtered operator

class FilteredC(BaseMatrix):
    def __init__(self, mata, matm_inv, tau, weights, freedofs = None):
        super().__init__()
        self.dt = tau
        self.weights = weights
        self.nsteps = len(weights)
        
        self.mata = mata
        self.matm_inv = matm_inv
        
        self.freedofs=freedofs
        
        
        self.vecu = self.mata.CreateColVector()
        self.tmpvec1 = self.mata.CreateColVector()
        self.tmpvec2 = self.mata.CreateColVector()


    def CreateColVector(self):
        return self.mata.CreateColVector()
    
    def Shape(self):
        return self.mata.shape
    
    def CreateVector(self,col):
        return self.mata.CreateVector(col)
    
    def Mult(self,rhs,out):
        with TaskManager():
            self.vecu.data = rhs
            tau = self.dt
            out.data = tau*weights[0]*self.vecu
            t = 0
            
            unew = self.tmpvec1
            uold = self.tmpvec2
            uold.data = self.vecu
            
            with TaskManager():
                for i in range(1,self.nsteps):
                    t += tau       
                    #print("\r time = {}, step = {}".format(t,i),end="")

                    unew.data = 2*self.vecu - uold 
                    unew.data -= tau**2 * self.matm_inv@self.mata * self.vecu
                    if self.freedofs:
                        unew.data[~self.freedofs] = 0.
                    uold.data = self.vecu
                    self.vecu.data = unew.data

                    out.data += tau*weights[i]*self.vecu

Geometry and mesh

r_l = 1.5
r_r = 0.15
d = 0.03
w = 0.03


circ_left = Circle( (-r_l-w,0), r_l).Face()
circ_right = Circle( (r_r+w,0), r_r).Face()
handle = MoveTo(-r_l,-d).Rectangle(r_l+r_r,2*d).Face()

dumbbell = circ_left+handle+circ_right
geo = OCCGeometry(dumbbell,dim=2)
Draw(dumbbell);

maxh = 2*d
mesh = Mesh(geo.GenerateMesh(maxh = maxh))
mesh.Curve(2)
Draw(mesh);

order = 2
fes = H1LumpingFESpace(mesh,order=order)
print("nDoF: ",fes.ndof)
u,v = fes.TnT()
mass = BilinearForm(u*v*dx(intrules=fes.GetIntegrationRules()),diagonal=True).Assemble().mat
stiffness = BilinearForm(grad(u)*grad(v)*dx).Assemble().mat
massinv = mass.Inverse()

nDoF:  14613

Determine time step by power iteration

tol = 1e-4
tmpv = stiffness.CreateVector()
tmpv2 = stiffness.CreateVector()
tmpv.SetRandom()
tmpv /= tmpv.Norm()
mus = [1]
for i in range(500):
    #print("i = {}".format(i))
    tmpv2.data = tmpv
    tmpv.data = massinv@stiffness*tmpv2
    mu = InnerProduct(tmpv2,tmpv)/InnerProduct(tmpv2,tmpv2)
    tmpv /= tmpv.Norm()
    mus.append(mu)
    if mus[-1]/mus[-2]<1+tol:
        break
        
lammax = mus[-1]*(1+tol*200)
print("estimated largest eigenvalue: ",lammax)
tau = np.sqrt(4/lammax)
print("by power iteration stable for tau = {}".format(tau))

estimated largest eigenvalue:  49228.18216580044
by power iteration stable for tau = 0.009014115036549765

Define and plot filter function

w_min = 0
w_max = 3

L = 200
endT = L*tau

weightf = lambda t: 4/np.pi*np.cos((w_max+w_min)/2*t)*(w_max-w_min)/2*np.sinc((w_max-w_min)/2*t/np.pi)

weights = weightf(tau*np.arange(L))


def beta(omega):
    if np.isscalar(omega):
        q = 1
    else:
        q = np.ones(omega.shape)
    
    q_old = q
    out = tau*weights[0]*q
    for alpha in weights[1:]:
        q_new = 2*q-tau**2*omega**2*q-q_old
        q_old = q
        q = q_new
        out += tau*alpha*q
    return out

omegas = np.arange(0,20,0.01)

betas = beta(omegas)
pl.plot(omegas,betas);

It remains to start the Krylov loop

maxsteps = 25


C = FilteredC(stiffness,massinv,tau,weights)


errsmin = []
resmin = []

tol = 1e-6
solveevery = 1

tmpvec = mass.CreateVector()
vecs = MultiVector(tmpvec,0)

tmpvec.SetRandom()
tmpvec.data = 1/tmpvec.Norm()*tmpvec

vecs.Append(tmpvec)


pl.figure()
pl.plot(omegas,betas);
pl.xlim((0,10))
for i in range(maxsteps):
    #print("\n Krylowstep = {}".format(i))

    tmpvec.data = C*vecs[-1]
    vecs.AppendOrthogonalize(tmpvec)
    
    if i%solveevery == 0:
        tvecs = MultiVector(tmpvec,len(vecs))
        tvecs.data = stiffness*vecs
        tvecs2 = MultiVector(tmpvec,len(vecs))
        tvecs2.data = mass*vecs

        matm_proj = InnerProduct(tvecs2,vecs)
        mats_proj = InnerProduct(tvecs,vecs)
        lam,v = spl.eigh(mats_proj.NumPy(),matm_proj.NumPy())

        lam = np.sqrt(np.abs(lam))
        for li in lam:
            pl.plot(li.real,1/maxsteps*i,'xr')
        

print("first 10 approximated eigenvalues")
for l in np.sort(lam)[:10]:
    print(l)

first 10 approximated eigenvalues
1.7068135620766916e-06
1.2015955703259844
1.2275073463471204
1.9051215442919904
2.0363109063162037
2.2840876468299025
2.579927270284838
2.801083497774604
2.889397348173861
3.5454948557875094

eigf = (vecs*Matrix(v.real)).Evaluate()

gfu = GridFunction(fes)

for i in range(1,5):
    gfu.vec.data = eigf[i]
    Draw(gfu)